function [ u_h condK ] = codinaSolveMatrix_BC_direct( bc_Type , N, id_bd, xNodes, s_nears, p_s, dp_s, xSamples, w_s )
%CODINASOLVEMATRIX_BC_DIRECT 
%
%   bc_Type  : The type of boundary condition to be used, 0 is homogeneus,
%              1 is standard Codina's problem
%   N        : Number of nodes - not used, keep for compatibility.
%   xNodes   : Node points.
%   xSamples : Sample points.
%   id_bd    : indices for the boundary nodes with dirichlet BC's.
%   w_s      : Gaussian integration weights for each sample point.
%   s_nears  : For each sample point, the nearest nodes.
%   p_s      : Shape function of each node evaluated at each sample point.
%   dp_s     : Shape function Gradient of each node evaluated at each
%              sample point.


% Number of node points
totalNodes      = size(xNodes,1);
% Number of sample points
totalSamples    = size(xSamples,1);
% Number of boundary nodes
nB = length(id_bd);

%% Stiffness matrix assembly-----------------------------------------------
K = zeros(totalNodes, totalNodes);
% For each sample point (Gauss Point)
for k = 1 : totalSamples;
    % Vector with the nearest node points of the k-ith gauss point.
    k_near = s_nears{k};
    % Gradient of the Shape function at the k-ith gauss point.
    dp_k   = dp_s{k};
    % Sum the contribution of the (\nabla u \cdot \nabla v)
    K(k_near,k_near) = K(k_near,k_near) + (dp_k*dp_k') * w_s(k);
end

%% Right hand side assembly------------------------------------------------
% Right hand side vector
rhs  = zeros(totalNodes, 1);

%% Homogeneus boundary condition problem
if (bc_Type == 0)
    % for each sample point k
    for k = 1 : totalSamples;
        % Gets the neares Nodes of the k-ith sample point.
        k_near = s_nears{k};
        % number of neares nodes of the k-ith sample point.
        n_k    = length(k_near);
        % Shape function value at the k-ith sample point.
        p_k    = p_s{k};
        % Gaussian integration weight for the k-ith sample point.
        w_gk   = w_s(k);
        % The external "force" (-2y)
        val    = -2.0 * xSamples(k,2) * w_gk;
        
        % loop on the neighbors of the i-th gauss points
        for ia = 1 : n_k;
            i      = k_near(ia);
            rhs(i) = rhs(i) + p_k(ia) * val;
        end;
    end;
    % impose the homogeneus boundary condition on the boundary
    for boundaryNodeIndex = 1 : nB;
        i      = id_bd(boundaryNodeIndex);
        rhs(i) = 0;
    end;
end;

%% Standard Codina's problem
if (bc_Type == 1);
  % Non-Homogeneous Boundary Conditions
  for boundaryNodeIndex = 1 : nB;
    nodeIndex  = id_bd(boundaryNodeIndex);
    % The Codina boundary condition: (y x (1 - x))
    b          = xNodes(nodeIndex,2) * ...
                 xNodes(nodeIndex,1) * (1 - xNodes(nodeIndex, 1));
    rhs(:)     = rhs(:) - b * K(:, nodeIndex);
  end;  
  
  for boundaryNodeIndex = 1 : nB;
    nodeIndex      = id_bd(boundaryNodeIndex);
    % The Codina boundary condition: (y x (1 - x))
    rhs(nodeIndex) = xNodes(nodeIndex,2) *...
                     xNodes(nodeIndex,1) * (1 - xNodes(nodeIndex,1));
  end;
end;

%% FEM lie imposition of boundary condition on the System matrix
% The rows and columns of each boundary node i are set to 0, except for the
% diagonal (i,i) which is set to 1.
for boundaryNodeIndex = 1 : nB
  nodeIndex              = id_bd(boundaryNodeIndex);
  K(nodeIndex,:)         = 0; 
  K(:,nodeIndex)         = 0;
  K(nodeIndex,nodeIndex) = 1;
end

%% Solve the linear system of equation.
u_h   = K \ rhs;
condK = cond(K);
end

